Learning vs Memory

Was the measured response interpreted as learning or memory?

mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_E1)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.3793   18.7586   28.7586   41.2577   29.4729   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0076  0.0873     30     no  Study_ID 
## sigma^2.2  0.0340  0.1843     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval   ci.lb 
## Learning_vs_memoryLearning    0.2303  0.0450  5.1227  90  <.0001  0.1410 
## Learning_vs_memoryMemory      0.1624  0.0355  4.5713  90  <.0001  0.0919 
##                              ci.ub 
## Learning_vs_memoryLearning  0.3197  *** 
## Learning_vs_memoryMemory    0.2330  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E1)

##   R2_marginal R2_coditional 
##    0.02572583    1.00000000

Figure corresponds to Fig. 4A in main text

LvsM_E <- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21) + # mean estimate
  scale_size_continuous(range = c(1, 7)) + # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

LvsM_E

Fig. S4a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E1 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = Learning_vs_memory,
    VCV = VCV_E)
res_table_mod_E1 <- get_estimate(model = mod_E1, contra = contra_mod_E1, moderator = Learning_vs_memory)

res_table_mod_E1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Learning          0.230    0.141    0.320  0.00000170   -0.185    0.645
## 2 Memory            0.162    0.0919   0.233  0.0000154    -0.249    0.574
## 3 Learning-Memory  -0.0679  -0.164    0.0284 0.165        -0.484    0.349

Type of assay

The broad category of the type of assay used to measure learning or memory.

dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)

mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)

summary(mod_E2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -7.8496   15.6993   27.6993   42.6311   28.7237   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0099  0.0995     30     no  Study_ID 
## sigma^2.2  0.0321  0.1792     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval    ci.lb   ci.ub 
## Type_assayConditioning    0.2167  0.0351  6.1650  89  <.0001   0.1468  0.2865 
## Type_assayHabituation     0.1821  0.1147  1.5886  89  0.1157  -0.0457  0.4100 
## Type_assayRecognition     0.0554  0.0640  0.8659  89  0.3889  -0.0717  0.1826 
##  
## Type_assayConditioning  *** 
## Type_assayHabituation 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E2)

##   R2_marginal R2_coditional 
##    0.07376134    1.00000000

Figure corresponds to Fig. 4B in main text

Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34", "grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Conditioning", "Recognition")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Learning_E

Fig. S4b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E2 <- contrast_fun(data = dat1, response = lnRR_Ea, moderator = Type_assay,
    VCV = VCV_E1)
res_table_mod_E2 <- get_estimate(model = mod_E2, contra = contra_mod_E2, moderator = Type_assay)

res_table_mod_E2

## # A tibble: 6 x 7
##   Levels                Estimate Lower_CI Upper_CI     P_value Lower_PI Upper_PI
##   <chr>                    <dbl>    <dbl>    <dbl>       <dbl>    <dbl>    <dbl>
## 1 Conditioning            0.217    0.147    0.287      2.02e-8   -0.197    0.630
## 2 Habituation             0.182   -0.0457   0.410      1.16e-1   -0.285    0.649
## 3 Recognition             0.0554  -0.0717   0.183      3.89e-1   -0.371    0.482
## 4 Conditioning-Habitua…  -0.0345  -0.262    0.193      7.63e-1   -0.501    0.432
## 5 Conditioning-Recogni…  -0.161   -0.296   -0.0265     1.96e-2   -0.590    0.268
## 6 Habituation-Recognit…  -0.127   -0.381    0.127      3.24e-1   -0.607    0.353

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~Type_reinforcement - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)

summary(mod_E3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -7.8702   15.7405   27.7405   42.6723   28.7649   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0109  0.1046     30     no  Study_ID 
## sigma^2.2  0.0319  0.1787     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
## 
## Model Results:
## 
##                                   estimate      se    tval  df    pval    ci.lb 
## Type_reinforcementAppetitive        0.2175  0.0726  2.9942  89  0.0036   0.0732 
## Type_reinforcementAversive          0.2191  0.0410  5.3456  89  <.0001   0.1377 
## Type_reinforcementNot applicable    0.0812  0.0560  1.4504  89  0.1505  -0.0300 
##                                    ci.ub 
## Type_reinforcementAppetitive      0.3618   ** 
## Type_reinforcementAversive        0.3006  *** 
## Type_reinforcementNot applicable  0.1924      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E3)

##   R2_marginal R2_coditional 
##    0.07720103    1.00000000

Figure corresponds to Fig. 4C in main text

Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34", "grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Aversive", "Not applicable")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Reinforcement_E

Fig. S4c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E3 <- contrast_fun(data = dat2, response = lnRR_Ea, moderator = Type_reinforcement,
    VCV = VCV_E2)
res_table_mod_E3 <- get_estimate(model = mod_E3, contra = contra_mod_E3, moderator = Type_reinforcement)

res_table_mod_E3

## # A tibble: 6 x 7
##   Levels                 Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>                     <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Appetitive              0.217     0.0732   0.362     3.56e-3   -0.218    0.653
## 2 Aversive                0.219     0.138    0.301     6.89e-7   -0.200    0.638
## 3 Not applicable          0.0812   -0.0300   0.192     1.50e-1   -0.345    0.507
## 4 Appetitive-Aversive     0.00164  -0.163    0.167     9.84e-1   -0.442    0.445
## 5 Appetitive-Not applic… -0.136    -0.315    0.0425    1.33e-1   -0.585    0.312
## 6 Aversive-Not applicab… -0.138    -0.259   -0.0172    2.56e-2   -0.567    0.291

Age at enrichment

The effect of age when individuals were exposed to environmental enrichment.

dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)


mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat3)

summary(mod_E4)

## 
## Multivariate Meta-Analysis Model (k = 88; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -6.9490   13.8980   23.8980   36.1698   24.6480   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0078  0.0885     29     no  Study_ID 
## sigma^2.2  0.0327  0.1808     88     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
## 
## Model Results:
## 
##                            estimate      se    tval  df    pval   ci.lb   ci.ub 
## Age_EE_exposureAdolescent    0.1799  0.0370  4.8553  86  <.0001  0.1062  0.2535 
## Age_EE_exposureAdult         0.2340  0.0620  3.7734  86  0.0003  0.1107  0.3573 
##  
## Age_EE_exposureAdolescent  *** 
## Age_EE_exposureAdult       *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E4)

##   R2_marginal R2_coditional 
##    0.01127347    1.00000000

Figure corresponds to Fig. 4D in main text

Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10))  

Age_E

Fig. S4d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E4 <- contrast_fun(data = dat3, response = lnRR_Ea, moderator = Age_EE_exposure,
    VCV = VCV_E3)
res_table_mod_E4 <- get_estimate(model = mod_E4, contra = contra_mod_E4, moderator = Age_EE_exposure)

res_table_mod_E4

## # A tibble: 3 x 7
##   Levels           Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>               <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Adolescent         0.180    0.106     0.254 0.00000533   -0.227    0.587
## 2 Adult              0.234    0.111     0.357 0.000295     -0.185    0.653
## 3 Adolescent-Adult   0.0541  -0.0895    0.198 0.456        -0.371    0.479

Exercise enrichment

Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?

mod_E5 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_E5)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -10.0303   20.0606   30.0606   42.5597   30.7749   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0096  0.0980     30     no  Study_ID 
## sigma^2.2  0.0345  0.1857     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_exerciseExercise       0.1849  0.0407  4.5464  90  <.0001  0.1041  0.2657 
## EE_exerciseNo exercise    0.1900  0.0556  3.4151  90  0.0010  0.0795  0.3005 
##  
## EE_exerciseExercise     *** 
## EE_exerciseNo exercise  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E5)

##   R2_marginal R2_coditional 
##  0.0001360993  0.9999999987

Figure corresponds to Fig. 4E in main text

Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Exercise_E

Fig. S4e Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E5 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = EE_exercise,
    VCV = VCV_E)
res_table_mod_E5 <- get_estimate(model = mod_E5, contra = contra_mod_E5, moderator = EE_exercise)

res_table_mod_E5

## # A tibble: 3 x 7
##   Levels               Estimate Lower_CI Upper_CI   P_value Lower_PI Upper_PI
##   <chr>                   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>    <dbl>
## 1 Exercise              0.185     0.104     0.266 0.0000169   -0.240    0.610
## 2 No exercise           0.190     0.0795    0.301 0.000958    -0.241    0.621
## 3 Exercise-No exercise  0.00511  -0.132     0.142 0.941       -0.434    0.444

Social enrichment

Did enrichment involve adding more conspecifics in the cage? (Note that we did not include studies that only provided social enrichment but no abiotic enrichment.)

dat4 <- filter(dat, EE_social %in% c("Social", "Non-social"))
VCV_E4 <- impute_covariance_matrix(vi = dat4$lnRRV_E, cluster = dat4$Study_ID, r = 0.5)

mod_E6 <- rma.mv(yi = lnRR_Ea, V = VCV_E4, mod = ~EE_social - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat4)

summary(mod_E6)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.3254   18.6507   28.6507   41.1498   29.3650   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0095  0.0974     30     no  Study_ID 
## sigma^2.2  0.0340  0.1844     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 801.5313, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 17.0110, p-val < .0001
## 
## Model Results:
## 
##                      estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_socialNon-social    0.1413  0.0498  2.8349  90  0.0057  0.0423  0.2403   ** 
## EE_socialSocial        0.2204  0.0432  5.0976  90  <.0001  0.1345  0.3063  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E6)

##   R2_marginal R2_coditional 
##    0.03377455    1.00000000

Figure corresponds to Fig. 4F in main text

Social_E <-orchard_plot(mod_E6, mod = "EE_social", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Social_E 

Fig. S4f Orchard plot showing the group-wise means of the categorical variable ‘EE_social’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_E6 <- contrast_fun(data = dat4, response = lnRR_Ea, moderator = EE_social,
    VCV = VCV_E4)
res_table_mod_E6 <- get_estimate(model = mod_E5, contra = contra_mod_E6, moderator = EE_social)

res_table_mod_E6

## # A tibble: 3 x 7
##   Levels                   Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 EE_exerciseExercise        0.185    0.104     0.266  1.69e-5   -0.240    0.610
## 2 EE_exerciseNo exercise     0.190    0.0795    0.301  9.58e-4   -0.241    0.621
## 3 EE_exerciseExercise-EE_…   0.0791  -0.0520    0.210  2.34e-1   -0.355    0.514

Multi-moderator model

We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.

# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
    filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
        c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
        "Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"))

# VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster =
# dat_Efm$Study_ID, r = 0.5)

mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay - 1 + Learning_vs_memory +
    Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Efm)
# summary(mod_Efm) r2_ml(mod_Efm)

res_Efm <- dredge(mod_Efm, trace = 2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))

The akaike weights for the top set of models with AIC < 6

dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2 <- subset(res_Efm, delta <= 6, recalc.weights = FALSE)
importance(res_Efm2)

##                      Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights:      0.63            0.36       0.16        0.14     
## N containing models:   11              10          5           5     
##                      Learning_vs_memory Type_reinforcement
## Sum of weights:      0.10               0.07              
## N containing models:    4                  4

Funnel plot

We used funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.

Figure corresponds to Fig. 7A in main text

# funnel plot
Funnel_E <- funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")

# Funnel_E

Fig.S5 Funnel plot showing the standard error and residuals (lnRR) from the full model.

Egger’s regression and time-lag bias

We fitted effective sample size calculated as:

\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]

and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.

VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID,
    r = 0.5)
dat_Efm$c_Year_published <- as.vector(scale(dat_Efm$Year_published, scale = F))

dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))

PB_MR_E <- rma.mv(lnRR_Sa, VCV_Efm, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory +
    Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise +
    Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
    test = "t", data = dat_Efm)

# PB_MR_E

estimates_PB_MR_E <- estimates.CI(PB_MR_E)
estimates_PB_MR_E

##                              estimate        mean        lower       upper
## 1                             intrcpt -9.48981044 -50.89863818 31.91901730
## 2                        sqrt_inv_e_n -0.06281042  -0.86857733  0.74295650
## 3            Learning_vs_memoryMemory -0.06858190  -0.17753636  0.04037255
## 4                      Year_published  0.00467568  -0.01591231  0.02526367
## 5               Type_assayHabituation -0.55467640  -0.99644706 -0.11290575
## 6               Type_assayRecognition -0.11960483  -0.45315288  0.21394322
## 7          Type_reinforcementAversive  0.15536913  -0.15356140  0.46429965
## 8    Type_reinforcementNot applicable  0.32741009  -0.11084714  0.76566731
## 9                     EE_socialSocial  0.04883160  -0.16642221  0.26408540
## 10             EE_exerciseNo exercise -0.03983198  -0.28097488  0.20131093
## 11           Age_stress_exposureAdult -0.16033127  -0.54043113  0.21976858
## 12 Age_stress_exposureEarly postnatal -0.05978991  -0.45195630  0.33237649
## 13        Age_stress_exposurePrenatal -0.13988907  -0.53943086  0.25965272

Leave-one-out analysis

We individually removed studies to determine if any singular studies have a disproportionate effect on the meta-analytic mean and CI.

dat$Study_ID <- as.factor(dat$Study_ID)

LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {

    d <- dat %>%
        filter(Study_ID != levels(dat$Study_ID)[i])

    VCV_Eb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)

    LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = VCV_Eb, random = list(~1 |
        Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
        levels(dat$Study_ID)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0) {
    df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
    return(df)
}


# using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
    bind_rows %>%
    mutate(left_out = levels(dat$Study_ID))


# saveRDS(MA_CVR_E,file = here('Rdata', 'MA_CVR_E.rds'))

Figure corresponds to Fig. S3 in Supplementary material S2.

# MA_CVR_E <- readRDS(file = here('Rdata', 'MA_CVR_E.rds'))

# telling ggplot to stop reordering factors
MA_CVR_E$left_out <- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out <- factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)


# plotting
leaveoneout_E <- ggplot(MA_CVR_E) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
    geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
    geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_E0$ci.ub,
    lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
    ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
    coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
    theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout_E

dat$Study_ID <- as.integer(dat$Study_ID)

Fig. S6 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.

Using SMD

We re-ran the main effect of environmental enrichment meta-analytic model using standardised mean difference (SMD) instead of lnRR.

mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), test = "t", data = dat)
summary(mod_E0a)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -111.6284   223.2568   231.2568   241.3003   231.7220   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.2512  0.5012     30     no  Study_ID 
## sigma^2.2  0.4247  0.6517     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.7241  0.1295  5.5895  91  <.0001  0.4668  0.9814  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

i2_ml(mod_E0a)

##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09

Risk of Bias

Testing a model including if the study was randomised as a moderator and conducting a pair-wise contrast analysis.

mod_randomised_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_randomisation - 1,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_randomised_E)
r2_ml(mod_randomised_E)

contra_mod_randomised_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_randomisation,
    VCV = VCV_E)
res_table_mod_randomised_E <- get_estimate(model = mod_randomised_E, contra = contra_mod_randomised_E,
    moderator = ROB_randomisation)

res_table_mod_randomised_E

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Unclear       0.140    0.0488    0.232 0.00304      -0.276    0.556
## 2 Yes           0.227    0.139     0.315 0.00000158   -0.188    0.642
## 3 Unclear-Yes   0.0868  -0.0352    0.209 0.161        -0.337    0.511

Testing a model including if the study used blinding as a moderator and conducting a pair-wise contrast analysis.

mod_blinding_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_blinding - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_blinding_E)
r2_ml(mod_blinding_E)

contra_mod_blinding_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_blinding,
    VCV = VCV_E)
res_table_mod_blinding_E <- get_estimate(model = mod_blinding_E, contra = contra_mod_blinding_E,
    moderator = ROB_blinding)

res_table_mod_blinding_E

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Unclear       0.198   0.122      0.274 0.00000131   -0.226    0.622
## 2 Yes           0.148   0.00811    0.288 0.0384       -0.292    0.588
## 3 Unclear-Yes  -0.0501 -0.208      0.108 0.530        -0.496    0.396

Learning vs Memory

Was the measured response interpreted as learning or memory?

mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_S1)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -14.5281   29.0562   39.0562   51.5553   39.7705   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0094  0.0970     30     no  Study_ID 
## sigma^2.2  0.0388  0.1969     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
## 
## Model Results:
## 
##                             estimate      se     tval  df    pval    ci.lb 
## Learning_vs_memoryLearning   -0.1211  0.0476  -2.5423  90  0.0127  -0.2157 
## Learning_vs_memoryMemory     -0.0974  0.0380  -2.5648  90  0.0120  -0.1728 
##                               ci.ub 
## Learning_vs_memoryLearning  -0.0265  * 
## Learning_vs_memoryMemory    -0.0219  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S1)

##   R2_marginal R2_coditional 
##   0.002776034   1.000000000

Figure corresponds to Fig. 5A in main text

LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

LvsM_S 

Fig. S8a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S1 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Learning_vs_memory,
    VCV = VCV_S)
res_table_mod_S1 <- get_estimate(model = mod_S1, contra = contra_mod_S1, moderator = Learning_vs_memory)

res_table_mod_S1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Learning         -0.121   -0.216   -0.0265  0.0127   -0.567    0.325
## 2 Memory           -0.0974  -0.173   -0.0219  0.0120   -0.540    0.345
## 3 Learning-Memory   0.0237  -0.0780   0.125   0.644    -0.424    0.472

Type of assay

The broad category of the type of assay used to measure learning or memory.

dat$Type_assay <- as.factor(dat$Type_assay)

VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)


mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)

summary(mod_S2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.8028   19.6055   31.6055   46.5374   32.6299   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0161  0.1271     30     no  Study_ID 
## sigma^2.2  0.0279  0.1671     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb    ci.ub 
## Type_assayConditioning   -0.0981  0.0375  -2.6192  89  0.0104  -0.1725  -0.0237 
## Type_assayHabituation    -0.4615  0.1126  -4.0969  89  <.0001  -0.6853  -0.2377 
## Type_assayRecognition    -0.0534  0.0645  -0.8287  89  0.4095  -0.1816   0.0747 
##  
## Type_assayConditioning    * 
## Type_assayHabituation   *** 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S2)

##   R2_marginal R2_coditional 
##     0.1853359     1.0000000

Figure corresponds to Fig. 5B in main text

Learning_S <- orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Habituation", "Recognition")),
              annotations = "***") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Learning_S

Fig. S8b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S2 <- contrast_fun(data = dat1, response = lnRR_Sa, moderator = Type_assay,
    VCV = VCV_S1)
res_table_mod_S2 <- get_estimate(model = mod_S2, contra = contra_mod_S2, moderator = Type_assay)

res_table_mod_S2

## # A tibble: 6 x 7
##   Levels                   Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Conditioning              -0.0981  -0.173   -0.0237  1.04e-2  -0.522    0.326 
## 2 Habituation               -0.461   -0.685   -0.238   9.21e-5  -0.935    0.0119
## 3 Recognition               -0.0534  -0.182    0.0747  4.10e-1  -0.490    0.383 
## 4 Conditioning-Habituation  -0.363   -0.584   -0.142   1.55e-3  -0.835    0.109 
## 5 Conditioning-Recognition   0.0447  -0.0882   0.177   5.06e-1  -0.393    0.482 
## 6 Habituation-Recognition    0.408    0.161    0.655   1.48e-3  -0.0768   0.893

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)

mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~Type_reinforcement - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)

summary(mod_S3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.8810   27.7621   39.7621   54.6939   40.7864   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0103  0.1016     30     no  Study_ID 
## sigma^2.2  0.0387  0.1966     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
## 
## Model Results:
## 
##                                   estimate      se     tval  df    pval 
## Type_reinforcementAppetitive       -0.1846  0.0749  -2.4649  89  0.0156 
## Type_reinforcementAversive         -0.0730  0.0427  -1.7081  89  0.0911 
## Type_reinforcementNot applicable   -0.1172  0.0590  -1.9851  89  0.0502 
##                                     ci.lb    ci.ub 
## Type_reinforcementAppetitive      -0.3334  -0.0358  * 
## Type_reinforcementAversive        -0.1579   0.0119  . 
## Type_reinforcementNot applicable  -0.2345   0.0001  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S3)

##   R2_marginal R2_coditional 
##     0.0366428     1.0000000

Figure corresponds to Fig. 5C in main text

Reinforcement_S <- orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Reinforcement_S

Fig. S8c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S3 <- contrast_fun(data = dat2, response = lnRR_Sa, moderator = Type_reinforcement,
    VCV = VCV_S2)
res_table_mod_S3 <- get_estimate(model = mod_S3, contra = contra_mod_S3, moderator = Type_reinforcement)

res_table_mod_S3

## # A tibble: 6 x 7
##   Levels                    Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                        <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Appetitive                 -0.185   -0.333  -3.58e-2  0.0156   -0.649    0.280
## 2 Aversive                   -0.0730  -0.158   1.19e-2  0.0911   -0.521    0.375
## 3 Not applicable             -0.117   -0.235   1.12e-4  0.0502   -0.572    0.338
## 4 Appetitive-Aversive         0.112   -0.0592  2.82e-1  0.198    -0.360    0.583
## 5 Appetitive-Not applicable   0.0674  -0.119   2.54e-1  0.474    -0.410    0.545
## 6 Aversive-Not applicable    -0.0442  -0.173   8.46e-2  0.497    -0.502    0.414

Age at stress

The age when individuals were exposed to stress.

mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_S4)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -12.4083   24.8166   38.8166   56.1579   40.2166   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0048  0.0694     30     no  Study_ID 
## sigma^2.2  0.0392  0.1979     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
## 
## Model Results:
## 
##                                     estimate      se     tval  df    pval 
## Age_stress_exposureAdolescent         0.0074  0.1159   0.0641  88  0.9490 
## Age_stress_exposureAdult             -0.2279  0.0622  -3.6664  88  0.0004 
## Age_stress_exposureEarly postnatal   -0.0561  0.0435  -1.2891  88  0.2007 
## Age_stress_exposurePrenatal          -0.1145  0.0743  -1.5404  88  0.1271 
##                                       ci.lb    ci.ub 
## Age_stress_exposureAdolescent       -0.2228   0.2377      
## Age_stress_exposureAdult            -0.3514  -0.1044  *** 
## Age_stress_exposureEarly postnatal  -0.1425   0.0304      
## Age_stress_exposurePrenatal         -0.2621   0.0332      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S4)

##   R2_marginal R2_coditional 
##    0.09987307    1.00000000

Figure corresponds to Fig. 5D in main text

Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
    geom_signif(comparisons = list(c("Adult", "Early postnatal")),
              annotations = "*") +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Age_S 

Fig. S8d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S4 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Age_stress_exposure,
    VCV = VCV_S)
res_table_mod_S4 <- get_estimate(model = mod_S4, contra = contra_mod_S4, moderator = Age_stress_exposure)

res_table_mod_S4

## # A tibble: 10 x 7
##    Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##    <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
##  1 Adolescent                0.00743  -0.223    0.238  9.49e-1   -0.469    0.484
##  2 Adult                    -0.228    -0.351   -0.104  4.20e-4   -0.663    0.207
##  3 Early postnatal          -0.0561   -0.142    0.0304 2.01e-1   -0.482    0.370
##  4 Prenatal                 -0.114    -0.262    0.0332 1.27e-1   -0.557    0.328
##  5 Adolescent-Adult         -0.235    -0.497    0.0260 7.70e-2   -0.727    0.257
##  6 Adolescent-Early postna… -0.0635   -0.309    0.182  6.09e-1   -0.547    0.421
##  7 Adolescent-Prenatal      -0.122    -0.395    0.152  3.78e-1   -0.620    0.377
##  8 Adult-Early postnatal     0.172     0.0211   0.323  2.60e-2   -0.271    0.615
##  9 Adult-Prenatal            0.113    -0.0791   0.306  2.45e-1   -0.346    0.573
## 10 Early postnatal-Prenatal -0.0584   -0.229    0.113  4.99e-1   -0.509    0.392

Type of stress

The type of stress used.

dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "Maternal separation",
    "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)

mod_S5 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat5)
summary(mod_S5)

## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -11.4138   22.8276   36.8276   53.5887   38.3618   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0115  0.1071     25     no  Study_ID 
## sigma^2.2  0.0401  0.2002     85     no     ES_ID 
## sigma^2.3  0.0000  0.0000      4     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
## 
## Model Results:
## 
##                                          estimate      se     tval  df    pval 
## Type_stress_exposureCombination           -0.0500  0.0892  -0.5605  81  0.5767 
## Type_stress_exposureMaternal separation   -0.0539  0.0560  -0.9630  81  0.3384 
## Type_stress_exposureNoise                 -0.1203  0.1036  -1.1608  81  0.2491 
## Type_stress_exposureRestraint             -0.2101  0.0704  -2.9863  81  0.0037 
##                                            ci.lb    ci.ub 
## Type_stress_exposureCombination          -0.2274   0.1274     
## Type_stress_exposureMaternal separation  -0.1653   0.0575     
## Type_stress_exposureNoise                -0.3265   0.0859     
## Type_stress_exposureRestraint            -0.3501  -0.0701  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S5)

##   R2_marginal R2_coditional 
##    0.07029554    1.00000000

Figure corresponds to Fig. 5E in main text

Stressor<- orchard_plot(mod_S5, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Stressor

Fig. S8e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S5 <- contrast_fun(data = dat5, response = lnRR_Sa, moderator = Type_stress_exposure,
    VCV = VCV_S3)
res_table_mod_S5 <- get_estimate(model = mod_S5, contra = contra_mod_S5, moderator = Type_stress_exposure)

res_table_mod_S5

## # A tibble: 10 x 7
##    Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##    <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
##  1 Combination              -0.0500    -0.227   0.127  0.577     -0.535    0.435
##  2 Maternal separation      -0.0539    -0.165   0.0575 0.338     -0.519    0.411
##  3 Noise                    -0.120     -0.326   0.0859 0.249     -0.617    0.376
##  4 Restraint                -0.210     -0.350  -0.0701 0.00373   -0.683    0.263
##  5 Combination-Maternal se… -0.00393   -0.213   0.206  0.970     -0.502    0.494
##  6 Combination-Noise        -0.0703    -0.342   0.202  0.608     -0.598    0.457
##  7 Combination-Restraint    -0.160     -0.386   0.0658 0.162     -0.665    0.345
##  8 Maternal separation-Noi… -0.0664    -0.301   0.168  0.575     -0.575    0.443
##  9 Maternal separation-Res… -0.156     -0.335   0.0227 0.0861    -0.642    0.330
## 10 Noise-Restraint          -0.0898    -0.339   0.159  0.475     -0.606    0.426

Stess duration

Was the stress acute or chronic?

dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)

mod_S6 <- rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat6)
summary(mod_S6)

## 
## Multivariate Meta-Analysis Model (k = 89; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.7898   27.5796   37.5796   49.9092   38.3204   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0062  0.0786     29     no  Study_ID 
## sigma^2.2  0.0391  0.1978     89     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb    ci.ub 
## Stress_durationAcute      0.0135  0.0659   0.2045  87  0.8384  -0.1176   0.1445 
## Stress_durationChronic   -0.1360  0.0373  -3.6456  87  0.0005  -0.2101  -0.0618 
##  
## Stress_durationAcute 
## Stress_durationChronic  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S6)

##   R2_marginal R2_coditional 
##    0.08245896    1.00000000

Figure corresponds to Fig. 5F in main text

Duration_S <- orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Duration_S 

Fig. S8f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_S6 <- contrast_fun(data = dat6, response = lnRR_Sa, moderator = Stress_duration,
    VCV = VCV_S4)
res_table_mod_S6 <- get_estimate(model = mod_S6, contra = contra_mod_S6, moderator = Stress_duration)

res_table_mod_S6

## # A tibble: 3 x 7
##   Levels        Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>            <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Acute           0.0135   -0.118  0.145   0.838      -0.429    0.456
## 2 Chronic        -0.136    -0.210 -0.0618  0.000453   -0.565    0.294
## 3 Acute-Chronic  -0.149    -0.300  0.00111 0.0517     -0.599    0.300

Multi-moderator model

We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.

# selecting moderator levels with k >=5
dat_Sfm <- dat %>%
    filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
        c("Appetitive", "Aversive", "Not applicable"), Type_stress_exposure %in%
        c("Restraint", "Noise", "Maternal separation", "Combination"), Stress_duration %in%
        c("Chronic", "Acute"))

# VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster =
# dat_Sfm$Study_ID, r = 0.5)

mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~Type_assay - 1 + Learning_vs_memory +
    Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Sfm)
# summary(mod_Sfm) r2_ml(mod_Sfm)

res_Sfm <- dredge(mod_Sfm, trace = 2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))

The akaike weights for the top set of models with AIC < 6

dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2 <- subset(res_Sfm, delta <= 6, recalc.weights = FALSE)
importance(res_Sfm2)

##                      Type_assay Stress_duration Type_stress_exposure
## Sum of weights:      0.91       0.88            0.22                
## N containing models:    7          6               2                
##                      Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights:      0.16               0.04                0.04              
## N containing models:    2                  1                   1

Funnel plot

We used funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.

Figure corresponds to Fig. 7B in main text.

# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")

# Funnel_S

Fig.S9 Funnel plot showing the standard error and residuals (lnRR) from the full model.

Egger’s regression and time-lag bias

We fitted effective sample size calculated as:

\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]

and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.

# calculating inv effective sample size for use in full meta-regression
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID,
    r = 0.5)

dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))

# time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))

PB_MR_S <- rma.mv(lnRR_Sa, VCV_Sfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
    Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
    Age_stress_exposure + Stress_duration, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), method = "REML", test = "t", data = dat_Sfm, control = list(optimizer = "optim",
    optmethod = "Nelder-Mead"))

# PB_MR_S

estimates_PB_MR_S <- estimates.CI(PB_MR_S)
estimates_PB_MR_S

##                                   estimate         mean       lower       upper
## 1                                  intrcpt  0.208757092 -0.47118467  0.88869885
## 2                             sqrt_inv_e_n  0.109877719 -0.90636737  1.12612281
## 3                         c_Year_published -0.006425821 -0.02804904  0.01519740
## 4                    Type_assayHabituation -0.663428076 -1.09349872 -0.23335743
## 5                    Type_assayRecognition -0.069789046 -0.41156145  0.27198335
## 6                 Learning_vs_memoryMemory -0.076712961 -0.19562432  0.04219839
## 7               Type_reinforcementAversive  0.039563987 -0.13761334  0.21674132
## 8         Type_reinforcementNot applicable  0.189481532 -0.16955283  0.54851589
## 9  Type_stress_exposureMaternal separation  0.233030567 -0.40478886  0.87084999
## 10               Type_stress_exposureNoise -0.020700600 -0.76667312  0.72527192
## 11           Type_stress_exposureRestraint -0.067490308 -0.43366600  0.29868539
## 12                Age_stress_exposureAdult -0.129875123 -0.51040393  0.25065368
## 13      Age_stress_exposureEarly postnatal -0.193461086 -0.89857848  0.51165631
## 14             Age_stress_exposurePrenatal -0.146406300 -0.57377820  0.28096560
## 15                  Stress_durationChronic -0.395711964 -0.65566557 -0.13575836

Leave-one-out sensitivity analysis

We individually removed studies to determine if any singular studies have a disproportionate effect on the meta-analytic mean and CI.

dat$Study_ID <- as.factor(dat$Study_ID)

LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
    d <- dat %>%
        filter(Study_ID != levels(dat$Study_ID)[i])

    VCV_Sb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)

    LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = VCV_Sb, random = list(~1 |
        Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
        levels(dat$Study_ID)[i], ])
}


# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_S0) {
    df <- data.frame(est = mod_S0$b, lower = mod_S0$ci.lb, upper = mod_S0$ci.ub)
    return(df)
}


# using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
    bind_rows %>%
    mutate(left_out = levels(dat$Study_ID))

# saveRDS(MA_CVR_S,file = here('Rdata', 'MA_CVR_S.rds'))

Figure corresponds to Fig. S4 in Supplementary material S2

# MA_CVR_S <- readRDS(file = here('Rdata', 'MA_CVR_S.rds'))

# telling ggplot to stop reordering factors
MA_CVR_S$left_out <- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out <- factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)

# plotting
leaveoneout_S <- ggplot(MA_CVR_S) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
    geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
    geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_S0$ci.ub,
    lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
    ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
    coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
    theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout_S

dat$Study_ID <- as.integer(dat$Study_ID)

Fig. S10 Leave-one-study-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.

Using SMD

We re-ran the main effect of stress meta-analytic model using standardised mean difference (SMD) instead of lnRR.

mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), test = "t", data = dat)
summary(mod_S0a)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -120.2817   240.5634   248.5634   258.6068   249.0285   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.3874  0.6224     30     no  Study_ID 
## sigma^2.2  0.4969  0.7049     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb    ci.ub 
##  -0.4256  0.1498  -2.8403  91  0.0056  -0.7232  -0.1279  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

i2_ml(mod_S0a)

##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09

Risk of Bias

Including if the study was randomised as a moderator and conducting a pair-wise contrast.

mod_randomised_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_randomisation - 1,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_randomised_S)
r2_ml(mod_randomised_S)

contra_mod_randomised_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_randomisation,
    VCV = VCV_S)
res_table_mod_randomised_S <- get_estimate(model = mod_randomised_S, contra = contra_mod_randomised_S,
    moderator = ROB_randomisation)

res_table_mod_randomised_S

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Unclear      -0.0672   -0.164   0.0291 0.169     -0.509    0.374
## 2 Yes          -0.139    -0.230  -0.0480 0.00317   -0.580    0.301
## 3 Unclear-Yes  -0.0719   -0.204   0.0607 0.284     -0.523    0.379

Including if the study used blinding and conducting a pair-wise contrast.

mod_blinding_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_blinding - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_blinding_S)
r2_ml(mod_blinding_S)

contra_mod_blinding_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_blinding,
    VCV = VCV_S)
res_table_mod_blinding_S <- get_estimate(model = mod_blinding_S, contra = contra_mod_blinding_S,
    moderator = ROB_blinding)

res_table_mod_blinding_S

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Unclear      -0.119    -0.196  -0.0419 0.00288   -0.564    0.326
## 2 Yes          -0.0567   -0.202   0.0890 0.442     -0.518    0.405
## 3 Unclear-Yes   0.0625   -0.102   0.227  0.454     -0.406    0.531

Learning vs Memory

Was the response interpreted as learning or memory?

mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_ES1)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -40.4769   80.9539   90.9539  103.4529   91.6682   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0292  0.1708     30     no  Study_ID 
## sigma^2.2  0.0232  0.1524     92     no     ES_ID 
## sigma^2.3  0.0034  0.0582      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval    ci.lb 
## Learning_vs_memoryLearning    0.1744  0.0722  2.4166  90  0.0177   0.0310 
## Learning_vs_memoryMemory      0.1057  0.0619  1.7065  90  0.0914  -0.0174 
##                              ci.ub 
## Learning_vs_memoryLearning  0.3178  * 
## Learning_vs_memoryMemory    0.2288  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES1)

##   R2_marginal R2_coditional 
##     0.0197648     0.9405080

Figure corresponds to Fig. 6A in main text

LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10))

LvsM_ES 

# axis.text.x = element_text(size = 10), # change font sizes
#         axis.text.y = element_text(size = 10),
#         legend.title = element_text(size = 7),
#         legend.text = element_text(size = 7)) 

Fig. S12a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES1 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Learning_vs_memory,
    VCV = VCV_ES)
res_table_mod_ES1 <- get_estimate(model = mod_ES1, contra = contra_mod_ES1, moderator = Learning_vs_memory)

res_table_mod_ES1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Learning          0.174    0.0310   0.318   0.0177   -0.316    0.665
## 2 Memory            0.106   -0.0174   0.229   0.0914   -0.379    0.591
## 3 Learning-Memory  -0.0687  -0.176    0.0382  0.205    -0.550    0.413

Type of assay

The broad category of the type of assay used to measure learning or memory.

VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)

mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)

summary(mod_ES2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0874   78.1748   90.1748  105.1066   91.1992   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0370  0.1924     30     no  Study_ID 
## sigma^2.2  0.0192  0.1386     92     no     ES_ID 
## sigma^2.3  0.0018  0.0422      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb   ci.ub 
## Type_assayConditioning    0.1525  0.0589   2.5877  89  0.0113   0.0354  0.2696 
## Type_assayHabituation     0.1990  0.1415   1.4070  89  0.1629  -0.0820  0.4801 
## Type_assayRecognition    -0.0048  0.0800  -0.0606  89  0.9518  -0.1637  0.1541 
##  
## Type_assayConditioning  * 
## Type_assayHabituation 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES2)

##   R2_marginal R2_coditional 
##    0.05775809    0.97105550

Figure corresponds to Fig. 6B in main text

Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Conditioning", "Recognition")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10)) 

Learning_ES

Fig. S12b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES2 <- contrast_fun(data = dat1, response = lnRR_ESa, moderator = Type_assay,
    VCV = VCV_ES1)
res_table_mod_ES2 <- get_estimate(model = mod_ES2, contra = contra_mod_ES2, moderator = Type_assay)

res_table_mod_ES2

## # A tibble: 6 x 7
##   Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Conditioning              0.153     0.0354   0.270   0.0113   -0.340    0.645
## 2 Habituation               0.199    -0.0820   0.480   0.163    -0.356    0.754
## 3 Recognition              -0.00485  -0.164    0.154   0.952    -0.509    0.499
## 4 Conditioning-Habituation  0.0465   -0.219    0.312   0.728    -0.501    0.594
## 5 Conditioning-Recognition -0.157    -0.304   -0.0105  0.0360   -0.658    0.343
## 6 Habituation-Recognition  -0.204    -0.488    0.0801  0.157    -0.760    0.353

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID,
    r = 0.5)

mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~Type_reinforcement - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)

summary(mod_ES3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0604   78.1208   90.1208  105.0526   91.1452   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0382  0.1954     30     no  Study_ID 
## sigma^2.2  0.0189  0.1377     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
## 
## Model Results:
## 
##                                   estimate      se    tval  df    pval    ci.lb 
## Type_reinforcementAppetitive        0.1007  0.1075  0.9366  89  0.3515  -0.1129 
## Type_reinforcementAversive          0.1573  0.0569  2.7618  89  0.0070   0.0441 
## Type_reinforcementNot applicable    0.0147  0.0702  0.2101  89  0.8341  -0.1247 
##                                    ci.ub 
## Type_reinforcementAppetitive      0.3143     
## Type_reinforcementAversive        0.2704  ** 
## Type_reinforcementNot applicable  0.1541     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES3)

##   R2_marginal R2_coditional 
##     0.0586952     1.0000000

Figure corresponds to Fig. 6C in main text

Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Aversive", "Not applicable")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10)) 

Reinforcement_ES 

Fig. S12c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES3 <- contrast_fun(data = dat2, response = lnRR_ESa, moderator = Type_reinforcement,
    VCV = VCV_ES2)
res_table_mod_ES3 <- get_estimate(model = mod_ES3, contra = contra_mod_ES3, moderator = Type_reinforcement)

res_table_mod_ES3

## # A tibble: 6 x 7
##   Levels                    Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                        <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Appetitive                  0.101   -0.113   0.314   0.351     -0.420    0.621
## 2 Aversive                    0.157    0.0441  0.270   0.00698   -0.331    0.645
## 3 Not applicable              0.0147  -0.125   0.154   0.834     -0.480    0.510
## 4 Appetitive-Aversive         0.0566  -0.182   0.295   0.639     -0.475    0.588
## 5 Appetitive-Not applicable  -0.0859  -0.333   0.161   0.492     -0.621    0.450
## 6 Aversive-Not applicable    -0.143   -0.278  -0.00744 0.0389    -0.636    0.351

Age at enrichment

The age when individuals were exposed to environmental enrichment.

VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID,
    r = 0.5)

mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat3)

summary(mod_ES4)

## 
## Multivariate Meta-Analysis Model (k = 88; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -36.6407   73.2813   83.2813   95.5531   84.0313   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0332  0.1822     29     no  Study_ID 
## sigma^2.2  0.0207  0.1439     88     no     ES_ID 
## sigma^2.3  0.0007  0.0265      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
## 
## Model Results:
## 
##                            estimate      se    tval  df    pval    ci.lb 
## Age_EE_exposureAdolescent    0.1284  0.0583  2.2006  86  0.0304   0.0124 
## Age_EE_exposureAdult         0.1247  0.0921  1.3538  86  0.1793  -0.0584 
##                             ci.ub 
## Age_EE_exposureAdolescent  0.2444  * 
## Age_EE_exposureAdult       0.3077    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES4)

##   R2_marginal R2_coditional 
##  0.0000400928  0.9871095822

Figure corresponds to Fig. 6D in main text

Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10)) 

Age_enrichment_ES

Fig. S12d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES4 <- contrast_fun(data = dat3, response = lnRR_ESa, moderator = Age_EE_exposure,
    VCV = VCV_ES3)
res_table_mod_ES4 <- get_estimate(model = mod_ES4, contra = contra_mod_ES4, moderator = Age_EE_exposure)

res_table_mod_ES4

## # A tibble: 3 x 7
##   Levels           Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>               <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Adolescent        0.128     0.0124    0.244  0.0304   -0.350    0.607
## 2 Adult             0.125    -0.0584    0.308  0.179    -0.375    0.624
## 3 Adolescent-Adult -0.00373  -0.213     0.205  0.972    -0.513    0.506

Age at stress

The age when individuals were exposed to stress.

mod_ES7 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_ES7)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0917   78.1834   92.1834  109.5247   93.5834   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0295  0.1718     30     no  Study_ID 
## sigma^2.2  0.0232  0.1522     92     no     ES_ID 
## sigma^2.3  0.0091  0.0954      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 286.4225, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.5932, p-val = 0.1832
## 
## Model Results:
## 
##                                     estimate      se     tval  df    pval 
## Age_stress_exposureAdolescent        -0.0137  0.1762  -0.0775  88  0.9384 
## Age_stress_exposureAdult              0.1677  0.1140   1.4708  88  0.1449 
## Age_stress_exposureEarly postnatal    0.1067  0.0920   1.1602  88  0.2491 
## Age_stress_exposurePrenatal           0.3179  0.1311   2.4254  88  0.0173 
##                                       ci.lb   ci.ub 
## Age_stress_exposureAdolescent       -0.3639  0.3366    
## Age_stress_exposureAdult            -0.0589  0.3942    
## Age_stress_exposureEarly postnatal  -0.0761  0.2896    
## Age_stress_exposurePrenatal          0.0574  0.5784  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES7)

##   R2_marginal R2_coditional 
##    0.09276574    0.86630830

Figure corresponds to Fig. 6E in main text

Age_stress_ES<-orchard_plot(mod_ES7, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10))

Age_stress_ES

Fig. S12e Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES7 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Age_stress_exposure,
    VCV = VCV_ES)
res_table_mod_ES7 <- get_estimate(model = mod_ES7, contra = contra_mod_ES7, moderator = Age_stress_exposure)

res_table_mod_ES7

## # A tibble: 10 x 7
##    Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##    <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
##  1 Adolescent                -0.0137  -0.364     0.337  0.938    -0.619    0.592
##  2 Adult                      0.168   -0.0589    0.394  0.145    -0.376    0.711
##  3 Early postnatal            0.107   -0.0761    0.290  0.249    -0.420    0.634
##  4 Prenatal                   0.318    0.0574    0.578  0.0173   -0.241    0.876
##  5 Adolescent-Adult           0.181   -0.236     0.598  0.390    -0.465    0.828
##  6 Adolescent-Early postna…   0.120   -0.275     0.515  0.546    -0.512    0.753
##  7 Adolescent-Prenatal        0.332   -0.105     0.768  0.135    -0.328    0.991
##  8 Adult-Early postnatal     -0.0609  -0.280     0.158  0.581    -0.601    0.479
##  9 Adult-Prenatal             0.150   -0.132     0.432  0.292    -0.419    0.719
## 10 Early postnatal-Prenatal   0.211   -0.0462    0.469  0.107    -0.346    0.768

Order to treatment exposure

The order in which individuals were exposed to enrichment and stress.

mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_ES10)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0408   78.0817   90.0817  105.0135   91.1061   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0316  0.1777     30     no  Study_ID 
## sigma^2.2  0.0227  0.1507     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 292.2561, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1546, p-val = 0.0287
## 
## Model Results:
## 
##                                 estimate      se     tval  df    pval    ci.lb 
## Exposure_orderConcurrently        0.1492  0.1208   1.2351  89  0.2201  -0.0909 
## Exposure_orderEnrichment first   -0.1782  0.1659  -1.0744  89  0.2856  -0.5079 
## Exposure_orderStress first        0.1370  0.0526   2.6046  89  0.0108   0.0325 
##                                  ci.ub 
## Exposure_orderConcurrently      0.3893    
## Exposure_orderEnrichment first  0.1514    
## Exposure_orderStress first      0.2414  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES10)

##   R2_marginal R2_coditional 
##     0.1027207     1.0000000

Figure corresponds to Fig. 6F in main text

Order_ES <- orchard_plot(mod_ES10, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 15), # change font sizes
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 10))

Order_ES 

Fig. S12f Orchard plot showing the group-wise means of the categorical variable ‘Exposure_order’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.

contra_mod_ES10 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Exposure_order,
    VCV = VCV_ES)
res_table_mod_ES10 <- get_estimate(model = mod_ES10, contra = contra_mod_ES10, moderator = Exposure_order)

res_table_mod_ES10